Layer III pyramidal cells in the prefrontal cortex reveal morphological changes in subjects with depression, schizophrenia, and suicide

Brodmann Area 46 (BA46) has long been regarded as a hotspot of disease pathology in individuals with schizophrenia (SCH) and major depressive disorder (MDD). Pyramidal neurons in layer III of the Brodmann Area 46 (BA46) project to other cortical regions and play a fundamental role in corticocortical and thalamocortical circuits. The AutoCUTS-LM pipeline was used to study the 3-dimensional structural morphology and spatial organization of pyramidal cells. Using quantitative light microscopy, we used stereology to calculate the entire volume of layer III in BA46 and the total number and density of pyramidal cells. Volume tensors estimated by the planar rotator quantified the volume, shape, and nucleus displacement of pyramidal cells. All of these assessments were carried out in four groups of subjects: controls (C, n = 10), SCH (n = 10), MDD (n = 8), and suicide subjects with a history of depression (SU, n = 11). SCH subjects had a significantly lower somal volume, total number, and density of pyramidal neurons when compared to C and tended to show a volume reduction in layer III of BA46. When comparing MDD subjects with C, the measured parameters were inclined to follow SCH, although there was only a significant reduction in pyramidal total cell number. While no morphometric differences were observed between SU and MDD, SU had a significantly higher total number of pyramidal cells and nucleus displacement than SCH. Finally, no differences in the spatial organization of pyramidal cells were found among groups. These results suggest that despite significant morphological alterations in layer III of BA46, which may impair prefrontal connections in people with SCH and MDD, the spatial organization of pyramidal cells remains the same across the four groups and suggests no defects in neuronal migration. The increased understanding of pyramidal cell biology may provide the cellular basis for symptoms and neuroimaging observations in SCH and MDD patients.

: BA46 was taken from four groups of subjects (Control (C), Major depressive disorder (MDD), Suicide (SU), and Schizophrenia (SCH)) who were matched for sex, age, hand, post mortem interval (PMI), and storage time. More information can be found in the paper by Chen et al. [1], where the same subjects were studied.

S2
Supplementary Note 2 3D-analysis -sample preparation The pyramidal cells in BA46 layer III were obtained from the extracted biopsies with a diameter of 3 mm. The direction of the biopsy should be such that the pial surface is perpendicular to the neocortex after sampling. Therefore, the neocortical layers were parallel to the cutting direction during tissue sampling and visible after staining. The entire process is depicted in Fig.1. The grey matter in BA46 is approximately 2.5 mm thick and covers all six neocortical layers. All tissues beyond 2 mm were removed from the sample using a razor knife since only layer III is of interest in this study. The biopsies were then sliced with a razor knife to a z-height of 1 mm and immersed in Phosphate-buffered saline (pH 7.3) with sucrose for one day before being rinsed twice at room temperature for 5 minutes each time in 0.05 Mol maleate buffer (pH 5.2). The samples were then processed and embedded in a Leica EM TP Automated Tissue Processor (Leica Microsystems, Brønshøj, Denmark). The biopsies were dehydrated using a graded ethanol sequence after one hour of staining with 1% uranyl acetate in maleate solution (70%, 86%, 96% and 99%, 20 min each). After dehydration, the sample was washed three times in ten minutes with 100% acetone and then infiltrated with 100% acetone/epon 1:1 overnight for 12 hours with continuous rotation. Infiltrated samples were incubated for one hour in pure Resin 812 before being placed in embedded molds and polymerized for 24 hours in a prewarmed oven (60 • ). The biopsies were positioned at the bottom of the embedded form with the pial surface perpendicular to the cutting direction of the knife. After the resin had cured completely, most of the white resin from the embedded sample was roughly trimmed using a high-speed milling machine (EM TRIM2, Leica) with an angle set at 60 • . The first section of the sample was stained with 1% toluidine and inspected under a microscope to detect layer III. Next, a glass knife was used to trim approximately a 1×1.4 mm 2 area, leaving only neurons in layers I-IV in the sample. Figure 1B portrays the tissue embedding procedure. The vertical axis (VA) was set to be perpendicular to the pial surface after the block of tissue containing BA46 was removed from the hemisphere. A paper ruler was inserted in the bottom of a box where the tissue block was mounted in 7% agarose. The block was sliced into uniform random 2.5 mm thick parallel vertical slabs S3 perpendicular to the anterior-posterior position using the SURS technique after the agarose had hardened. The tissue block was then reassembled, and every second the slab was systematically sampled. Anatomical landmarks and VA of each slab were then characterized and documented, including the superior frontal gyrus, medial frontal gyrus, and inferior frontal gyrus. The slabs were then numbered and placed in a tissue processing embedding cassette with foam pads to keep them safe while submerged in a solution. The tissue was then dehydrated in a Shandon Citadel 1000 tissue processor (Thermo Scientific, USA) using a 70%, 96%, and 99% ethanol series. Next, the tissue was submerged in glycol methacrylate for two shifts in 1:1 ethanol: Technovit 7100 solution, then in pure Technovit 7100 without the hardener for one shift. Each shift lasted one day in a 4 • fridge. The tissue was then placed into an embedding mold after being withdrawn from the tissue processing embedding cassette. The tissue was put at the bottom of each embedding mold, which contained 15 ml of Technovit 7100. After that, 1 mL of hardener was applied, and the tissue took about five days to polymerize fully. The slabs were then cut into 40 µm sections, stained with 1% Toluidinblue-Borax Solution for 30 minutes, submerged in distilled water for 2 minutes, air-dried, mounted with Eukitt, and covered with 120 µm thick coverslips.

Stereological analysis
Every section was scanned using an Olympus scanning microscope with an Olympus 20x oil lens prior to stereological examination. The image viewer application OlyVIA v.3.2 (Olympus, USA) can display a high-resolution image with the option to zoom in and out of the slice, making it simpler to detect and deliniate BA46 based on the surrounding cytoarchitecture. The stereological study was then carried out for each subject in layer III of the whole BA46 using an Olympus BX51 light microscope with an Olympus DP70 camera (Olympus, Tokyo, Japan), an H138A motorized stage (Prior, Cambridge, UK) to load the microscopy slides, and stereological software VIS version 2017.7.1.3832 (Visiopharm, Hørsholm, Denmark). The amount of field of view (FOV) for each section was decided by the proportion of tissue that required to be examined to produce an accurate approximation, see table S2. This percentage was calculated by doing a cost-effective analysis based on pilot counts, with the objective of assessing an average of 2-300 points per participant. Layer III of BA46 was, therefore, ready to be sampled and examined using stereological probe sets in a systematic manner. Layer III of BA46 was then sampled and analyzed using stereological probes(Cavalieri estimator, optical fractionator probe, planar rotator). All S4 histological sections were counted while the investigator (NYL) was blinded to avoid bias.  Figure S1: Biopsies were used for 3D analysis, while the remaining portion of the tissue containing BA46 was used for stereological analysis.

S6
A B Figure S2: (A) A 7x7 point grid (green cross) was superimposed on the top of the section with a magnification of 1.25x. Layer III was within the region of interest (ROI) and was marked within green dashed lines. Points within the ROI were counted and represented with red circles. Scale bar= 2000 µm. (B) A field of view (FOV) with an x60 lens, which was superimposed on a counting frame (CF) on the top of the section. Only those cells that were within the box or touched the green line were counted. Cells that were outside of the box or hit the red line were excluded. Scale bar= 35 µm.

B A
iii ii iv i Figure S3: (A) Steps were used to sample a pyramidal cell with the Section Estimator. Within the optical disector, a neuron in focus (i). The blue line represents the predefined vertical axis, with the nucleolus serving as the reference point (ii). The top and bottom cell boundaries were marked (iii). There are four half-lines perpendicular to the vertical axis that appear in a uniform random position (iv). The + sign denotes the point where the neuron border and the half-line cross. Scale bar= 35 µm. (B) 2D representation of the displacement vectorc from the nucleolus (reference point) to the center of mass and the Miles ellipsoidē.ē represents the population shape. A red circle represents a particle's center of mass, while a black circle represents the reference point.

Supplementary Note 3 Point pattern analysis
The data consists of 39 point patterns containing 3D coordinates for the locations of centroids of pyramidal cells in layer III of BA46. Such a point pattern is considered as a realisation of a stochastic mechanism called a point process. Each point pattern corresponds to a subject, and each subject belongs to one of the four groups, control, depression, schizophrenia, and suicide. There are 10 subjects in the control group (C), 8 subjects who suffered from major depressive disorder (MDD), 10 subjects who suffered from schizophrenia (SCH), and 11 subjects who had committed suicide with a history of depression (SU).

The cylindrical K-function
To detect columnar structures in the point patterns, we used the cylindrical K-function, which was introduced in [2], and has previously been used to analyze the spatial structure of pyramidal cells (see e.g. [3], [4],[5]). The cylindrical K-function depends on a direction u, a radius r, and a height t. We use the notation K u (r, t) for the cylindrical K-function and ρ for the intensity of the underlying point process (the expected number of points per unit volume, see Fig.S5). To understand the intuitive interpretation of K u (r, t), consider a cylinder C v in 3D space with direction u, base radius r, and height 2t which is centred at a randomly selected point v of the underlying point process X. Then, ρK u (r, t) is the expected number of further points from X within C v . If a point process exhibits a columnar structure in a specific direction, K u (r, t) is expected to be particularly high for some range of r and t values when u is that direction. For estimating K u (r, t) from a point pattern, we use the nonparametric estimate from [2]. For the choice of u, we considered the directions of the three main axes; for the choice of r, we considered 128 equidistant points in the interval from 0µm to 25µm; and for t we used t = 80. Note that when trying to detect a columnar structure, it makes sense to use r < t.
An important assumption for using the cylindrical K-function is that the underlying point process is homogeneous. We assessed this assumption from the observed point patterns by looking at the histograms of each coordinate and 2D projections. All point patterns appeared reasonably homogeneous.
To use the estimatesK u (r, t) of the cylindrical K-function to make statements about columnar structures in the data, we compared the obtained estimates with the situation of a so-called homogeneous Poisson process also referred to as 'complete spatial randomness' (CSR). As the name suggests, S16 CSR is the situation where there is no structure in data. We make the comparison based on a number of simulated point patterns under CSR. The estimate of the cylindrical K-function is calculated for each of these simulations and the resulting curves are ordered by means of the so-called extreme rank length measure [6,7]. The ordered curves can then be used to construct a global envelope test for the null hypothesis that the observed point pattern is a realisation of CSR as described in [6]. We used a 95% global envelope test and 2000 simulations under CSR (the number of simulations followed the recommendations in the above references). The method returns a p-value and a graphical interpretation in the form of a global envelope which consists of a region in which the observed curve falls completely if and only if the global envelope test cannot be rejected at approximately level 5%. If the empirical curve falls above the envelope for some range of r-values, it indicates that the estimate is higher than expected under CSR, which in turn indicates some clustering. If this behaviour is most pronounced for a specific direction u, it indicates cylinder-shaped clusters in this direction. If the curve falls below the envelope, it is smaller than expected under CSR, which indicates the repulsion between points. Figures S6, S7, S8, and S9 show the cylindrical K-function for all three directions for each subject plus 95% global envelopes under CSR and the p-values of the corresponding tests. For the majority of subjects, and in particular for all 10 control subjects, 7 in MDD, 7 in SCH, and 3 in SU, there is clear evidence of a columnar structure in the direction of the x-axis, especially when considering radii between 5 and 20µm. In some of the point patterns, there is some evidence of repulsive behavior between points, which is very intuitive since each cell occupies an amount of physical space, which is not accounted for, and the cells cannot overlap. However, this may not explain all repulsion. To get a better visualisation for comparison, Figure S10 shows all empirical cylindrical K-functions in the same plot split by direction and group. Notice that there is a lot of variation in the empirical cylindrical K-function between subjects, and it is difficult to spot whether there is a clear difference between groups. To summarise the general behaviour of the empirical cylindrical K-function within each group, we pool the estimates in the following way. Let K ij u (r, t) be the empirical cylindrical K-function of the i ′ th point pattern in the j ′ th group and let this point pattern have n ij points. Let there be m j point patterns in group j. To estimate the cylindrical K-function for group j, we use the weighted mean This weighted mean was recommended for estimating Ripley's K-function (see below) from replicated point patterns in [8] (both when the windows are the same and different sizes and when there is a common and pattern-specific S22 intensity), and the argumentation seems reasonable for the cylindrical Kfunction as well. Figure S11 shows these estimates for each group in the three considered directions. This confirms that there is an overall tendency for the cylindrical K-function in the direction of the x-axis to reach particularly high values, which indicates a columnar structure in this direction. To decide whether the apparent differences between the empirical cylindrical K-functions estimated for each group are significant, we use random permutations to construct a global envelope test to get an approximate test of the null hypothesis that all point patterns are realisations of spatial point processes with the same cylindrical K-function in the following way. Denote the number of groups by g. By concatenatingK 1 u , . . . ,K g u , one curve of data is achieved. This curve is then compared to the situation under the null hy-S23 pothesis by making a number of random permutations of the point patterns; recalculatingK j u (r, t) for each group, which no longer consists of the same point patterns but has the same number of point patterns; concatenating the estimates; and finally using these curves to construct a global envelope test as above. We considered a 95% global envelope test based on 8000 random permutations. The observed concatenated curve will fall completely within the envelope if and only if the null hypothesis of no difference in the cylindrical K-function between groups cannot be rejected at level 5%. If the concatenated curves fall outside the envelope for any r-value, the test is rejected, and there are significant differences between the groups. The results of such tests for the cylindrical K-function directed along one of the three main axes can be seen in Figures S12, S13, and S14 where we for better visualisation split the concatenated curves into one plot for each group. These tests demonstrate that the cylindrical K-function along the x, y, and z-axes did not detect any significant difference between the groups.

Ripley's K-function
Ripley's K-function is a popular functional summary statistic for spatial point processes. It depends on an argument r, and we denote it by K(r). If ρ is the intensity of the process, ρK(r) is interpreted as the expected number of further points within the distance r of a typical point of the point process. Note that this interpretation is similar to the interpretation of the cylindrical K-function, where we just consider a cylinder instead of a sphere. Ripley's Kfunction may detect some different structures in data than the cylindrical Kfunction. If K(r) is higher (lower) than it is under CSR, it suggests that the point process is more clustered (repulsive) at interpoint distances r. When S27 using Ripley's K-function, we again assume homogeneity. For estimating K(r) from the data, we used the translation correction estimator [9, 10] and we compared these estimates to the case of CSR with global envelopes as described for the cylindrical K-function . Figures S15, S16, S17, and S18 show the estimates minus the theoretical value under CSR and 95% global envelope for each subject. In all point patterns we see some deviations from CSR. Most point patterns show repulsive behaviour for a range of r-values, which is again very intuitive because the cells do not overlap. However, since we also see a repulsion for radii much higher than the approximate diameter of a cell, which is about 10µm, it suggests that not all repulsion can be explained by the fact that the cells do not overlap. We also see evidence of clustering in about three-quarters of the subjects, especially for radii in the range of 5 − 10µm.   To summarise the empirical curves from each group, we used (1) with the estimate of the cylindrical K-function replaced by the estimate of Ripley's K-function, and we thus follow the recommendation in [8]. Figure S20 shows the resulting curves, which look very similar across groups. To test whether there are significant differences in Ripley's K-function between the groups, we used a 95% global envelope test in the same way as ex-S32 plained for the cylindrical K-function. The result can be seen in Figure S21, which shows no significant differences in Ripley's K-function between groups.

The empty-space function F
The empty-space function is a functional summary statistic of spatial point processes which is based on interpoint distances. It depends on an argument r, and we use the notation F (r). The function F (r) is interpreted as the probability of finding a point from the process within a sphere of radius r centered at a fixed location in space. When using the F -function, we again assume homogeneity. For estimating F (r) from the data, we used the reduced sample estimator [9] and we compared these estimates to the case S33 of CSR with global envelopes as described for the cylindrical K-function. Figures S22, S23, S24, and S25 show the estimates and 95% global envelopes for each subject. The envelopes are in all cases very narrow. The F -function detects deviations from CSR in 28 of the point patterns.    To get an estimateF j (r) of the F -function for the j'th group, we make the weighted average in equation (19) of [9]. It only makes sense to consider one estimate of F for each group if all point patterns within a group are assumed to be realisations of point processes with the same intensity. Based on Figure S5, this appears to be a reasonable assumption. Note that this assumption was not necessary in the corresponding analyses using the cylindrical K-function and Ripley's K-function. Figure S27 shows the resulting curves, which do not show much difference between groups. To test whether there are significant differences in the F -function between the groups, we again used a 95% global envelope test in the same way as explained for the cylindrical K-function. The result can be seen in Figure S28, which shows no significant differences.

The nearest-neighbour function G
The nearest-neighbour function is again a functional summary statistic of spatial point processes which is based on interpoint distances. It depends on an argument r, and we use the notation G(r). Intuitively, G(r) can be interpreted as the probability of finding another point of the process within a sphere of radius r centered at a random point of the process. Thus, if G(r) is higher (lower) than it is under CSR, it suggests that the point process is more clustered (repulsive) at interpoint distances r. When using the Gfunction, we again assume homogeneity. For estimating G(r) from the data, we used the reduced sample estimator [9] and we compared these estimates to the case of CSR with global envelopes as described for the cylindrical K-     To get an estimateĜ j (r) of the G-function for the j'th group, we make the weighted averageĜ ij as recommended in [11]. HereĜ ij (r) is the estimate of the G-function obtained from the i'th point pattern in the j'th group, and the remaining notation is as in (1). It again only makes sense to consider one estimate of G for each group if all point patterns within a group are assumed to be realisations of point processes with the same intensity. Figure S34 shows the resulting curves, which do not show much difference between groups. To test whether there are significant differences in the G-function between the groups, we again used a 95% global envelope test in the same way as explained for the cylindrical K-function. The result can be seen in Figure S35, which shows no significant differences. In spatial point pattern analyses, it is also very common to consider the functional summary statistic J(r) = (1−G(r))/(1−F (r)), which is constantly equal to 1 in the case of CSR. This very simple expression under CSR makes it very easy to make visual comparisons between an estimate of J and the theoretical value under CSR. However, we do not consider this summary statistic for the following reason. The estimation method used to estimate F for a 3D point pattern gives an approximation to F (r) which should not be compared to its theoretical value under CSR [see 9]. Since the estimate of F does not behave as its theoretical value under CSR, neither will the estimate of J and we thus lose the easy visual comparison with CSR. Therefore, we do not consider J.